The number of significant digits may be imprecise, prefer to use relative error.
In [ ]:
# Definition 1: Round down to p-sig. digit number
x = 0.90
x1 = 0.99 # 2 correct significant digit, actual difference 0.09
x2 = 0.89 # 1 correct significant digit, actual difference 0.01
In [ ]:
# Definition 2: Round to nearest p-sig. digit number
y = 0.9951 # --> 0.10
y1 = 0.9499 # --> 0.90 , only 1 correct sig. digit
y2 = 1.0000 # --> 0.10 , 3 correct sig. digits
The last definition avoids issues that the previous two had, but it's still non-ideal
In [ ]:
# Definition 3: Numbers x and x' match to p-sig. digits if x - x' < half a unit in p-th sig. digit of x
x1 = 0.123
x2 = 0.127
# 0.004 < (0.01 / 2) => x1 and x2 match in 2 significant digits according to this definition wchich may be slightly confusing
In [ ]:
def absolute_error(true_value, approx_value):
return abs(true_value - approx_value)
print 'Absolute error: {0:.9f}'.format(absolute_error(10.951, 10.949))
def relative_error(true_value, approx_value):
return absolute_error(true_value, approx_value) / abs(true_value)
print 'Relative error: {0:.9f}'.format(relative_error(10.951, 10.949))
In [ ]:
import numpy as np
def relative_error(true_value, approx_value):
return np.max(np.fabs((true_value - approx_value) / true_value))
x_value = np.array([10000, 0.01])
x_approx = np.array([9999, 0.02])
print relative_error(x_value, x_approx)
An important example of an ill-conditioned problem is finding roots of a polynomial. Let's look at Wilkinson’s polynomial as an example.
In [ ]:
# Wilkinson's polynomial is defined as p(x) = (x - 1)(x - 2)...(x - 20)
import numpy as np
x = np.linspace(0, 20, 4000)
y = 1
for i in range(1, 20):
y *= (x - i)
import matplotlib.pyplot as plt
%matplotlib inline
plt.ylim([-4e11, 4e11])
plt.plot(x, y)
Quoting Wikipedia:
Error coming from representing a function or continuous variable using finite number of evaluations - outside of scope of this notebook, mentioned for completeness only.
In [ ]:
import numpy as np
x = np.float64(0.1)
y = np.float32(0.1)
print '32- vs 64-bit representation difference:', abs(x - y)
In [ ]:
import numpy as np
# For explanation why we we use np.round rather than default Python 2.7.3 round function, see below
x = 9.945309
print np.round(x, 2), np.round(np.round(x, 2), 1)
print np.round(x, 1)
There are few differences between built-in Python 2.7 round function and numpy (a)round:
Note that Python 3 has a different round function that behaves more similarly to numpy round.
In [ ]:
import numpy as np
for i in range(13):
x = -3 + 0.5 * i
print '\t{0:5.1f}\t{1:5.1f}\t{2:5.1f}'.format(x, round(x), np.round(x))
In [ ]:
x_value = 0.123123
y_value = 0.123000
# We want to learn the value of d = x - y:
d_value = x_value - y_value
print 'Actual d value:',
# Assuming we're apprximating above calculation with a 4 decimal digits precision:
import numpy as np
d_approx = np.round(x_value, 4) - np.round(y_value, 4)
print 'Approx d value:', d_approx
print 'Absolute error: {0:.9f}'.format(abs(d_value - d_approx))
print 'Relative error: {0:.9f}'.format(abs((d_value - d_approx) / d_value))
Cancellation is an example of loss of significance and it happens when two nearly equal numbers are subtracted and can lead to significant inaccuracies. As an example let's look at function:
f(x) = (1 - cos x) / x^2
We can see that calculating f(x) near zero may lead to issues, since cos(0) = 1.
In [ ]:
import numpy as np
from math import sin, cos
near_zero = 1.2e-8
def f(x):
return (1 - np.cos(x)) / x**2
print 'Value of f near zero:', f(near_zero)
x = np.linspace(near_zero, 10, 100)
y = f(x)
import matplotlib.pyplot as plt
%matplotlib inline
plt.ylim([0, 1])
plt.plot(x, y)
Fortunately, we can rewrite f(x) in a form that is less prone to cancellation:
In [ ]:
def g(x):
return 0.5 * (2 * np.sin(x / 2) / x)**2
print 'Value of g near zero:', g(near_zero)
x = np.linspace(near_zero, 10, 100)
y = g(x)
import matplotlib.pyplot as plt
%matplotlib inline
plt.ylim([0, 1])
plt.plot(x, y)
Summing numbers: https://docs.python.org/2/library/math.html https://en.m.wikipedia.org/wiki/Kahan_summation_algorithm Interestingly, IPython Notebook appears to do the right thing by default, what? IPython session:
In [6]: from math import fsum
In [7]: sum([0.1] * 10) Out[7]: 0.9999999999999999
In [8]: fsum([0.1] * 10) Out[8]: 1.0
In [ ]:
from math import fsum
print '{0:0.20f}'.format(sum([0.1] * 10))
print '{0:0.20f}'.format(fsum([0.1] * 10))
In [ ]:
def naive_e(n):
return (1 + 1.0 / n)**n
from math import exp
e = exp(1)
for i in range(5, 20):
print naive_e(10**i) - e
In [ ]:
from math import sqrt
def identity(x, n):
for i in xrange(n):
x = sqrt(x)
for i in xrange(n):
x = x**2
return x
x = 2
for i in xrange(35, 60):
print x - identity(x, i)
In [ ]:
# Computing f(x) = (exp(x) - 1) / x == sum(x^i / (i + 1)!)
from math import exp, log
def f1(x):
if 0 == x:
return 1
return (exp(x) - 1) / x
def f2(x):
if 0 == x:
return 1
y = exp(x)
return (y - 1) / log(y)
# f(epsilon) ~= 1
for i in range(8, 15):
epsilon = 1.0 / (10**i)
print 'epsilon:', epsilon
print '|1 - f1(epsilon)|:', abs(1 - f1(epsilon))
print '|1 - f2(epsilon)|:', abs(1 - f2(epsilon))
print
# NOTE: Above doesn't hold if we calculate for powers of 2!
# for i in range(30, 40):
# print abs(1 - f1(1.0 / (2**i))), abs(1 - f2(1.0 / (2**i)))
In [ ]:
import numpy as np
def r(x):
# Calculate value of a rational function using Horner's rule (https://en.wikipedia.org/wiki/Horner%27s_method)
# More on the function can be looked up on Wolfram's Alpha:
# http://www.wolframalpha.com/input/?i=f(x)+%3D+(622.0+-+x+*+(751.0+-+x+*+(324.0+-+x+*+(59.0+-+4+*+x))))+%2F+(112+-+x+*+(151+-+x+*+(72+-+x+*+(14+-+x))))&t=crmtb01
p = 622.0 - x * (751.0 - x * (324.0 - x * (59.0 - 4 * x)))
q = 112 - x * (151 - x * (72 - x * (14 - x)))
return p / q
def calc_f(a):
t = np.array([a + k * 2**-52 for k in xrange(400)])
t = r(t)
t -= t[0]
t *= 1.0 / max(abs(t))
return t
import matplotlib.pyplot as plt
%matplotlib inline
def plot(t):
plt.plot(t, linestyle='--', marker='o')
plt.show()
for a in [1.606, 4, 8, 16, 32]:
plot(calc_f(a))
There are two concepts that we refer to as variance:
Variance of a random variable X is defined as:
Var(X) = E[(X - E[X])^2] = ... = E[X^2] - (E[X])^2
Where E[Z] is the expected value of a random variable Z.
The second form (E[X^2] - (E[X])^2) should be avoided when performing calculations on a fixed precision machine. Although it has a nice property that for sample variance it is easy to implement it while traversing the data just once, naive implementations usually suffer from extreme cancellation.
In [ ]:
# Calculate the variance of X = sum of dice in N throws. The probabilities of each side of dice are given.
import numpy as np
# Dice sides probabilities, p[i] = probability of throwing i. We use fractions to be able to compare our calculations to
# 'exact' values:
from fractions import Fraction as F
from numpy.random import randint as ri
max = 1000000
pf = np.array([0] + [F(ri(max), 1) for _ in range(6)])
pf /= sum(pf)
p = np.array([float(f) for f in pf])
# Number of throws
N = 10000 # 30000
# Dynamic program, we're holding probability of throwing k-dice for each possibility.
# First iteration:
dp = np.ones(1)
for i in range(N):
dp = np.convolve(dp, p)
print dp
# Let's calculate variance using both ways.
dice = np.arange(len(dp))
# var1 = E[(X - E[X])^2]
ex = (p * np.arange(7)).sum() * N
var1 = (dp * (dice - ex)**2).sum()
print 'Variance calculated using: E[(X - E[X])^2]: ', var1
# var2 = E[X^2] - (E[X])^2
ex2 = (dice**2 * dp).sum()
var2 = ex2 - ex**2
print 'Variance calculated using: E[X^2] - (E[X])^2: ', var2
print 'E[X^2]: ', ex2
print '(E[X2])^2: ', ex**2
# There is a simpler way to calculate variance in this particular problem, since all N throws are independent we can simply
# calculate variability of one and multiply it by N. To make sure this calculation is as precise as possible, we use Python
# representation of rational numbers.
# Let's use equation 1: E[(X - E[X])^2]
rangef = np.array([F(i, 1) for i in range(7)])
exf = (rangef * pf).sum()
ex2f = (rangef**2 * pf).sum()
varf = ex2f - exf**2
varf *= N
print 'Variance calculated fractions: ', float(varf)
print 'Relative error of: E[(X - E[X])^2]: ', abs((var1 - varf) / varf)
print 'Relative error of: E[X^2] - (E[X])^2: ', abs((var2 - varf) / varf)
There's more to calculating variance, but a good rule of a thumb is - if you're fine calculating variance based on E[(X - E[X])^2] equation, you should be fine.